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We study the effect of a static electric field on lithium donor spins in silicon. The anisotropy of 
the effective mass leads to the anisotropy of the quadratic Stark susceptibility, which we determined 
using the Dalgarno- Lewis exact summation method. The theory is asymptotically exact in the field 
domain below Li-donor ionization threshold, relevant to the Stark-tuning electron spin resonance 
experiments. To obtain the generalized Stark susceptibilities at arbitrary fields, we propose a new 
variational wave function which reproduces the exact results in the low-field limit. With the calcu- 
lated susceptibilities at hand, we are able to predict and analyze several important physical effects. 
First, we observe that the energy level shifts due to the quadratic Stark effect for Li donors in 
Si are equivalent to, and can be mapped onto, those produced by an external stress. Second, we 
demonstrate that the Stark effect anisotropy, combined with the unique valley-orbit splitting of a 
Li donor in Si, spin-orbit interaction and specially tuned external stress, may lead to a very strong 
modulation of the donor spin p-factor by the electric field. Third, we investigate the influence of 
random strains on the g-factor shifts and quantify the random strain limits and requirements to Si 
material purity necessary to observe the <?-factor Stark shifts experimentally. Finally, we discuss 
possible implications of our results for quantum information processing with Li spin qubits in Si. 
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I. INTRODUCTION 

Recent research has shown that electron spins bound to 
donors in Si are some of the most viable candidates for 
scalable, solid-state quantum computing applications 1 . 
The lithium donor is unique among other Si shallow 
donors^ because its ground state is degenerated We 
will demonstrate that the inverted electronic structure 
of the lithium donor allows intriguing possibilities for 
the efficient manipulation of Li spin qubits using stress 
and electric fields. The ability to couple individual elec- 
tron spins with local electric fields is a key feature of 
many spin-control proposals for quantum computing and 
single-electron spintronic applications J^r— The electri- 
cal control of single electron spins based on various cou- 
pling mechanisms has been demonstrated experimentally 
for semiconductor quantum dotsJ^r— Also, a possibil- 
ity of spin-orbit-driven electrical modulation of shallow 
donor g-factors in GaAs has been predicted theoreti- 
cally^ In this work, we present theoretical evidence that 
the Li donor in Si is yet another shallow donor system 
that possesses such functionality. This system may play 
a significant role in silicon-based quantum information 
processing (QIP) with spin qubits. 

The Stark effect for shallow donors in silicon has been 
studied both experimentally^ and theoretically^ - — 
Most of the theoretical calculations have been performed 
for substitutional donors (e.g. P)2ii but not for Li donors 
occupying tetrahedral interstitials in a silicon latticed 
Furthermore, most of the previous studies of substitu- 
tional donors ignore the interplay between the Stark and 
Zeeman effects. This interplay occurs because the spin- 
orbit interaction couples the non-degenerate ground state 
of a substitutional donor to the rest of Is manifok h 23 ' 24 



which is separated form the ground state by large valley- 
orbit splitting, A vo ^10 meV. Since the spin-orbit cou- 
pling is much smaller than A VOl the electrical modulation 
of the ^-factors of substitutional donors is strongly sup- 
pressed. On the contrary, the ground state of the inter- 
stitial Li donor is degenerate and even modest spin-orbit 
interaction may have a profound effect on the energy lev- 
els leading to a strong electric field dependence of the 
Zeeman splittings. 

Intriguingly, the orbital degeneracy of the ground state 
quintet of the lithium donor gives rise not only to non- 
trivial spin-orbit effects^ but also to a strong long-range 
elastic- dipole coupling between the orbital states of dif- 
ferent lithium donors^ The elastic-dipole coupling was 
studied earlier in connection to an acceptor-based quan- 
tum computing scheme in silicon proposed by Dykman 
and Golding 2 ^. Thus, it is of interest to study the lithium 
donor electron as a new candidate for QIP applications 
with a long-range inter-donor coupling. 

Experimental techniques utilizing pulsed electron spin 
resonance (ESR) measurements in interdigitated devices 
based on Sb-doped Si revealed observable Stark shifts 
of ESR lines in an electric field £. The shifts can 
be parametrized as AA/A = r] a £ 2 and Ag/g = r] g £ 2 , 
for the hyperfine constant and the ^-factor shifts re- 
spectively, where r\ a w —3.7 x 10 -3 ^tm 2 /V 2 and r/ g w 
— 1 x 10 -5 /im 2 /V 2 f^ We anticipate that experiments on 
similar devices based on Li-doped Si may reveal substan- 
tially larger shifts of the magnetic resonance lines on the 
order of 10 gauss for electric fields ranging from to 
3 kV/cm. This proposition seems somewhat surprising 
as there is no hyperfine interaction in the ground state 
for Li donors and the strength of the spin-orbit inter- 
action on the Li atom in Si is negligible^ Nonetheless, 
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spin-orbit effects have been observed quite prominently 
in ESR spectra of Li-doped natural Si under external 
stress^ 

This happens because the unique electronic structure 
of a Li donor enables observation of spin-orbit effects 
which have crystalline rather than impurity origin^ The 
situation in P-doped Si is quite different because the 
crystalline spin-orbit interaction is much less important 
due to its weak influence on the isolated, non-degenerate 
ground state. More specifically, the orbital states of the 
Is manifold of a shallow donor in silicon form a singlet 
A±, a doublet E and a triplet T2. For substitutional 
donors the singlet ground state A\ is separated by a gap 
exceeding lOmeV from the closely spaced doublet and 
triplet. The sequence of levels is inverted for the inter- 
stitial Li donors in such a way that the ground state 
is five-fold degenerate and composed of the orbital dou- 
blet ls(E) and triplet ls(T2) while the singlet ls(Ai) lies 
1.76 meV above (see Fig. Q]). 
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FIG. 1. Upper table shows the energy level diagrams for 
Is electron of the shallow donors in silicon: left column 
corresponds to EMT calculations giving the binding energy 
31.27 meV^ Experimental values for the Is energies of sub- 
stitutional phosphorus^ and interstitial lithium^ donors are 
shown in right and middle columns, respectively. Ionization 
energy of lithium donor 33.02 meV is in particularly close 
agreement with EMT result. Inset depicts the locations of 
6 conduction band minima in silicon. Lower table shows the 
experimental results for the energy differences between the 
singlet and triplet (A at) and between the doublet and triplet 
(Aet)- 

Within the Effective Mass Theory (EMT) the Is donor 
electron manifold has 6-fold orbital degeneracy corre- 
sponding to 6 minima of the silicon conduction band k SJ - . 
These minima are located along coordinate axes at about 
85% of the distance to the zone boundary as shown in the 
inset of Fig. [TJ 

h S j=sKQ rij, s = ±1, ko ~ 0.85 Gioo/2, (1) 

where rij is a Cartesian unit vector, j = x,y or z, Gooi = 
47r/ao is a magnitude of the reciprocal lattice vector in 
[001] direction, and ao is the Si lattice constant. The 



position of each valley is characterized by a composite 
index {sj}, where s = ±1 describes the valley centered on 
either the positive or the negative semi-axis, respectively. 

A short-range tetrahedral potential corresponding to 
the local symmetry of the donor site splits the valley 
degeneracy. The intervalley effects^ are described by 
the valley-orbit Hamiltonian 

H vo = {E + A ) \ si ) H + A i H I s *) (~ s * I 

z,s s,i 

+ A 2 (l-SiMWjl (2) 

i, j 7 8,s' 

Here, Eq is the binding energy corresponding to the so- 
lution of the single- valley Coulomb problem and the pa- 
rameters Ai are the matrix elements of the short-range 
central cell potential in the basis of the six valley orbitals 

(r\sj) = cxp (ik sj ■ r) u sj (r)F s:j (r), (3) 

where u S j (r) is a periodic part of the Bloch function cor- 
responding to the center of the valley k SJ , and F S j(r) is 
the Is envelope function which, along with Eq, is found 
by solving the single- valley Coulomb problem. 

The six eigenstates of the Hamiltonian ([2]) , referred to 
as symmetrized valley orbitals, can be expressed as 

lM>=E<l s J>' M= 6. (4) 

Each orbital in Eq. (0| belongs to the irreducible rep- 
resentation /i of the tetrahedral group Td characterized 
by the valley-orbit coefficients ct^ that are given in Ta- 
ble n 

The singlet state A\ and the doublet states Eg, E e 
are "even" with respect to the axis inversion, i.e. they 
are symmetric combinations of the opposite valley or- 
bitals, a S j=a- S j. On the contrary, each triplet state 
T^j (j=x,y or z) is an "odd" antisymmetric combination 
of just two opposite valleys orbitals, a S j=— a- S j. We 
have labeled the symmetrized valley orbitals according to 
Watkins and Ham^ and indicated their transformational 
properties under the group Td (last column in Table Uj). 
Thus the singlet, triplet and doublet states possess s, p 
and d-like characters, respectively. 

For phosphorus donors in silicon, Friesen2£ discussed 
two competing mechanisms which determine the behav- 
ior of the electron levels in the Is manifold. First, as 
the electric field is increased, the quadratic Stark shift 
causes the energy of each valley to decrease. This will 
influence the form of the valley-orbit Hamiltonian 
Formally, we have to replace Eq with valley-dependent 
terms E^ = Eq — 5i(£) and place these under the sum- 
mation sign in Eq. $2$, where the corrections Si include 
the quadratic Stark shifts of each valley. Secondly, simi- 
lar substitutions are required for the central cell terms in 
the Hamiltonian Ao — > Aoi, Ai — > An and A2 — > Azij. 
These parameters A are the matrix elements of the short- 
range central cell potential between the states \is); they 
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TABLE I. Table of valley-orbit coefficients for silicon donor Is 
orbital states. First column defines the order of states. Sec- 
ond column shows the state label derived from the correspond- 
ing representation of the group Td- Each row in the third 
column gives a set of 6 coefficients for sj=±l,±2,±3 

corresponding to the silicon conduction valley minima at 
±x, ±y, ±z axes, respectively. Fourth column shows basis 
functions describing the transformational properties of the 
corresponding states under the operations of the group Td 



are proportional to the amplitudes of the envelope func- 
tions F s j(0) in the central cell. Since the electric field £ 
pulls the electron away from the donor site, the ampli- 
tudes F S j(0) and the matrix elements A are decreasing 
functions of £ 2 . As a result, the energy spectrum of the 
manifold narrows and the ground state shifts upward in 
the direction opposite to the quadratic Stark shift. 

Our studies reveal that the narrowing effect is not as 
important for Li as for P donors. Because the valley-orbit 
splitting for Li is considerably smaller than phosphorus, 
the narrowing of the spectrum does not overwhelm the 
quadratic Stark effect in determining the overall behav- 
ior of the Is manifold in the presence of an electric field. 
According to our findings, the most important effect for 
Li is the anisotropy of the quadratic Stark effect. This 
anisotropy allows the electric field to induce unique split- 
ting of the Li ground state and leads to a very non-trivial 
interplay of the Zccman and Stark effects. 

Our approach utilizes the Dalgarno-Lewis exact sum- 
mation method to determine the quadratic Stark suscep- 
tibility. While this calculation is quite lengthy, it pro- 
vides us with an important calibration tool to further de- 
vise our variational function in the presence of the electric 
field and gauge it against the exact small-field asymptotic 
behavior. We will discuss how the ground state splitting 
caused by the electric field can replicate stress effects and 
can potentially be used to manipulate a Li spin qubit. 
While most of the theoretical papers on the Stark effect 
were concerned with rather large fields ; 21 ' 29 " — our stud- 
ies are focused on the field domain below the 3-5kV/cm 
relevant to the ESR Stark experiments, i.e. below ion- 
ization threshold of shallow donors in Si. 

The paper is organized as follows. In Section |TT] 
we study the quadratic Stark effect for a single-valley 
Schrodinger equation for a shallow donor in Si. The 
goal of this section is to calculate the Stark suscepti- 
bility, which is asymptotically exact in the limit of low 
electric fields using an exact summation method of Dal- 



garno and Lewis tailored to account for the effective mass 
anisotropy. Based on the findings of Section [IIJ we pro- 
pose a new variational wave function in Section [IIII This 
wave function not only replicates the exact susceptibility 
at low fields but also describes the off-center displacement 
of the probability maximum for intermediate and high 
electric fields. In Sections llVltVl we introduce valley-orbit 
and external stress effects, and in Section fVTl we describe 
the spin-orbit and Zeeman Hamiltonians. In Section [VIII 
based on the results of the previous sections, we calculate 
the electric-field-induced ESR g-factor shifts for various 
types of spin-flip transitions and analyze the effects of 
random strains on the Stark shifts of the ESR spectra. 
Section [Villi contains the summary and conclusions. 



II. QUADRATIC STARK EFFECT 

As a first step let us consider a quadratic Stark ef- 
fect for a single- valley donor. We start with the EMT 
Kohn-Luttingcr Hamiltonian--i perturbed by an external 
electric field £ 



H- 



n 2 



d 2 



f d 2 

2mi \dx 2 dy 2 ^ dz 2 



e 



e£r, (5) 



where k is the dielectric constant, 7 = m^/m^ is the ef- 
fective mass anisotropy parameter, m± and my are trans- 
verse and longitudinal effective masses respectively, and 
we assume that the heavy-mass axis of the valley is along 
z. For our purposes it is convenient to rewrite the Hamil- 
tonian ([5]) using a scaling transformation z — > yFyz. This 
yields 



H- 



if- 



2m± \dx 2 



d 2 . 

dy 2 



dz 2 



e 

h'p 



e£ p, (6) 



where p = (x,y,^z). 

Our immediate goal is to find a second order shift of 
the ground state energy produced by an external static 
electric field £ (Stark shift): 



1 



(ls\e£ ■ p\n)(n\e£ ■ p\ls) 



Eis — E n 



— -^Xal3£ a £/3, 



(7) 



where we introduced the Stark susceptibility tensor and 
assumed summation over repeating Greek indices. In the 
chosen coordinate system, the tensor \ a p is diagonal and 
has only two distinct components Xzz = Xll an d Xxx = 
Xyy = X-L due to the axial symmetry of the valley with 
respect to z. 

It is well known from the classical problem of the 
quadratic Stark effect in a hydrogen atom that any per- 
turbative treatment of the hydrogen Schrodinger equa- 
tion in an external electric field requires summation of 
infinite series to account for the excited states of the con- 



tinuous spectrum 
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To accomplish this for the lithium 
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donor, we will utilize the Dalgarno-Lewis exact summa- 
tion method 3 - 6 .. We define a vector function /(r) to sat- 
sify the equation 



[f(r),H}\ls) = -p\ls) 



(8) 



With this equation in hand, consider an important con- 
dition for computing the second-order summation: 



(n\p\ls) = -(n\[f,H]\ls) 

= {E n -E ls )(n\f\ls) 



(9) 



As a consequence, if the solution /(r) of Eq. ([8]) is known, 
the Stark susceptibility tensor defined in Eq. ([7]) can be 
obtained as 



X a p = 2e 2 (ls\r a ■ /a(r)|ls>. 



(10) 



As both, the potential terms of H, as well as /(r) 
are functions of the coordinates only, they will commute 
and thus [f,H] = -(ti 2 /2m_ L ) [/, V 2 ] . Expressed in the 
coordinate-representation, Eq. © becomes a differential 
equation defining the function /(r) such that 



V 2 /(r) + 2(|- V)/(r) 



2mj_ 



where £ = Vln(V'is), and ipi s {r) = (r| Is) is the ground 
state wave function. 

To describe the ground state wave function tpi s {r) 
we use the Kohn-Luttinger functio n 26 ' 28 with ellipsoidal 
symmetry, 



^ ls ( r ) = , 3 .1/2 ex P I _ : I ' ( 12 ) 



(Trai) 



where /3 = ja^/a 2 and a_L,ay are the transverse and 
longitudinal radii of the isolated valley in the absence of 
the electric field. Then 



1 xx + yy + fizz 
a_L ^Jx 2 + y 2 + f3z 2 



(13) 



At this point it is convenient to introduce another scal- 
ing transformation such that x — > a±x, y — > a±y, and 
z — > aj_z/^/]3. Under this transformation all the coor- 
dinates become dimcnsionless (measured in units of a±) 
and the ground state wave function ipi s — > exp(— r)/ ^fir. 



Also, in Eq. ([TT]) we explicitly separate the anisotropic 
term proportional to A = 1 — /3 and treat it as a per- 
turbation (i.e. we will construct our solution as a power 
series in the anisotropy parameter A). This yields: 



(14) 



where 



and 



D r = V 2 - 28/dr, (15a) 
t) z = d 2 /dz 2 - 2(z/r)d/dz. (15b) 
From the axial symmetry of Eq. (|14|) . we are looking 



for solutions in the form: 



f x (r) = {a±/E ± )ooa<f>-f±(r,0), 
f y (r) = (a x /E 1 _)smcj>-f x (r,9), 
f z (r) = (a ]] /E ± )-f ll (r,6), 



(16a) 
(16b) 
(16c) 



where Ej_ = h /2m±a^_ and ay = aj_y7//3. The par- 
tial differential equation (|14p may be further simplified 
by generating a coupled system of ordinary differential 
equations. The explicit angular dependence of the oper- 
ators D r and D z is presented and analyzed in Appendix 
|A1 The analysis suggests that /||(r, 9) and /j_(r, 9) can be 
expanded into series in Legendre's and associated Legen- 
dre's polynomials, respectively, as follows: 



f\\{r,9)=Y,f\\,i{r)Pi{cos9), 
i=i 

/ ± (r,0) = -^/ ± ,(r)F i 1 (cos( 



(17a) 
(17b) 



i=i 



Note that the series expansion begins at I = 1. 

Substituting the Legendre's expansions (|17a[) and 
(|17b| into Eqs (|16a|) - (|16c|) . inserting the results into 
Eq. (|10p for the Stark susceptibility tensor, and integrat- 
ing over the angular variables we obtain: 



XII 



3£_L Jo 

8e 2 a 2 L f 00 
3E ± J a 



r 3 exp(-2r)/| !a (r)dr, (18a) 
r 3 exp(-2r)/ ± ,i(r)dr, (18b) 



The final task is to compute the radial parts f±^(r) and 
f\\ t l(r) of our function /. We accomplish this by inserting 
the combination of Eqs (|16a[) - (|16c[) and (|17a|) - (|17b|) into 
Eq. (|14[) . This procedure results in the following system 
of ordinary differential equations for the radial functions 
defined when / > (see Appendix [S] for more details) : 



J 



5 



r 2 D r f u - 1(1 + 1)/, M = A(d?/| M+2 + /3°/, M +7°/||, i - 2 ) - r 3 S n , 
r 2 D r f±.i - 1(1 + l)f±,i = Ha]f±,i+2 + $lf±,i + lif±,l-2) - r 3 S ni 

where a™ , /3 ; m , 7™ are second order radial differential operators having the common structure 

d 2 



(0) 



dr 2 



Hi 



9 , (2) , (3) 
- 2r I T -Q r + a l,m r + a l,m 



(19a) 
(19b) 



(20) 



with different coefficients oifln, /3; and 7^ as specified 
in Appendix El Note that 7J — 7 J = 0. 

The solution to Eq. (fT9|) with A = is known (see Ap- 
pendix This suggests that a perturbative treatment 
of the differential problem is appropriate. With this in 
mind, we expand f± i(r) and /n as power series in 



I 

A = 1 — p. Letting q = {|| , _L}, our radial functions take 
the form: 



4Kr) = /S ) W + A/ 9 7(r) + A^;/(r) + 



(21) 



Equating terms of like order in A, we obtain a chain of 
differential equations for different I > 1 and n > 0: 



r 2 D r f, 



(«) 



1(1 + = (i _ «5„ o) (af/fc? + A /||" _1) + Tf/fcl?) - ^ 3 5n<5„o 



T>f> r f$ i(i + = (i - m te/<y + 2 + AViT 15 + 7^/17-2 



a J-(n-l) 



(22a) 
(22b) 



From this point forward, the solution method will be 
entirely iterative. Determination of requires f^, 

fqi requires f^l and so forth. Note that Eqs. (f22|) with 
n = reduce to Eqs. (fT"9|) with A = 0. Since we know 
the solution for A = (sec Appendix [A]) . we use it as 

zeroth iteration, solve each differential equation for /^y 
and proceed to the higher orders. The first few solutions 



.(0) _ ,(0) _ r(r + 2) 
J\\,l - J±,i 4 

{1) _ r(r + 3) 

J \\;l 5 ' 

(1) _ r(r - 2) 
J±A - 40 ■ 



(23a) 
(23b) 
(23c) 



the expanded /, 
8e 2 a 2 



X± = 



3E ± 



8e 2 aj 
3E X 



(0.84375 + 0.8250A + 0.8173A 2 + ...), 



(24a) 



(0.84375 + 0.0094A + 0.0058A 2 + ...), 



(24b) 



where for silicon mi = .191m e , mn = .916m e , a±_ = 
23.65 A, and an = 13.60 A. For these values of the 
effective masses and Bohr radii E± = 35.68 meV, A = 
0.370, and our susceptibilities are 



X|| « 1.74/zeV(kV/cm)- 2 
X± « 3.54^eV(kV/cm)- 2 . 

III. VARIATIONAL METHOD 



(25a) 
(25b) 



The analytic form of the solutions past this point quickly 
becomes cumbersome. However, graphing the integrands 
A™r 3 exp(— 2r)fg " (r) of the first three terms shows the 
rapid convergence of the series, as seen in Fig. 

The final expressions for the susceptibilities, utilizing 



The findings of Scction[IT]provide very accurate asymp- 
totic behavior of the energy levels at low electric fields. 
To achieve high accuracy for the Stark shifts of the ESR 
lines both at low and intermediate fields ~5 kV/cm we 
have to extend the scope of our methodology and use 
a specially crafted variational approach that is valid at 
higher fields and replicates the low-field results of the 
previous section. 
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FIG. 2. (Color online.) The figures show successive inte- 
grands of (fT8l) up to the second order term for fields lying in 
the ^-direction (a) and the z-direction (b). The vertical axes 
are scaled appropriately for each graph, and the legend for 
both figures is shown in (b). The inset in (a) is a magnified 
view of the two dashed lines (n = 1,2) that are not resolved 
in the main plot. The inset in (b) displays dipole moments 
induced by the electric fields parallel and perpendicular to the 
heavy-mass axis z. 



Thus we seek to devise and use a trial function which 
will adequately represent the perturbed donor state, ac- 
curate up to at least the second order. This variational 
function will be used to extract expansions of the en- 
ergy and the central cell electron density, which cannot 
be obtained from the second order perturbation theory. 
At low fields we will see converging results between the 
variational method of this section and the infinite-series 
perturbation theory of Section [TTJ 

In EMT, we can represent the orbital states of donors 



in silicon with hydrogcnic-like envelope functions. We 
expect that a homogeneous electric field will admix p- 
state components into this ground state. In hydrogen, 
each p-function with principle quantum number n will 
be of the form rL^ l _ 2 (r/n)exp(—2r/n) where E\_ 2 is 
the generalized Laguerrc polynomial of degree n — 2. The 
higher p-states will contribute terms of order r n ~ l to the 
perturbed state. 

The higher p-states with principle quantum numbers 
n > 3 couple somewhat weaker to the ground state. How- 
ever, since there are many of these states, it is reasonable 
to surmise that any highly accurate approximate method 
must account for their contributions. In summary, to ac- 
curately reflect the influence of the electric field on the 
hydrogenic wave functions of donors, our trial function 
must take into account the higher order radial contribu- 
tions from the excited states. 

It is instructive to illustrate our approach by consider- 
ing the standard Hamiltonian for a hydrogen atom placed 
in a homogeneous electric field parallel to z-axis: 



H 



2m„ 



e 
r 



e£z. 



(26) 



It can be shown that the exact first order expansion of the 
wave function, and the associated second order expansion 
of the ground state energy are 



,32,33 



1 p(r) = Mr)+e£f(r)Mr), (27) 
E=-E Ry -e 2 £ 2 (^ \zf(r)\^ Q ) 1 (28) 

where f(r) = (l/aBE Ry )(r + 2as)z/4 with as and E Ry 
arc the Bohr radius and the Rydberg energy respec- 
tively. The function f(r) = fj(r) with j = x,y, z, is the 
isotropic analog of the function f(r) considered in the 
infinite-series summation of the previous section. The 
wave function (l27l) can be recast as 



(29) 



ip(r) = [1 + (qi + q 2 r) z] exp(-r/a), 



where qi = e£/2E Ry , q 2 = e£ j 4aBE Ry , and a — as- 

To extend our formalism to higher orders we can iden- 
tify qi and q 2 as the first order expansion of some un- 
known variational parameters. Therefore we assume that 
the trial function takes the form of Eq. (|29|) with the un- 
known variational parameters qi = qi(£), <?2 = <72 (£)) 
and a = a{£) to be determined via a standard mini- 
mization routine at an arbitrary field £. This procedure 
leads to the following expansion of the energy expecta- 
tion value: 



(E) = 



= -Ef 



9 



£ 2 /E R y 



(30) 



which is in complete agreement with the exact second- 
order £-field expansion of the ground state energy (|28|) . 

The term proportional to rz in Eq. (|29[) is responsible 
for admixture of the excited states into the ground state 
by the electric field. By including this term in the trial 
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function we were able to replicate the exact second or- 
der correction to the energy. In addition, the variational 
method allows for an efficient (albeit approximate) ac- 
count of higher order terms outside of the perturbative 
regime. For the hydrogen atom, setting q 2 = in Eq. (|29|) 
yields the following Taylor expansion of the expectation 
value of energy: 

(E) = -E Ry - a 2 B e 2 £ 2 /E Ry + ... . (31) 

Thus by neglecting contributions of the excited states, 
we would introduce a relative error of 11% to the second 
order energy shift, running the risk of producing more 
significant errors in the higher orders. 

In the case of the silicon donor states, our varia- 
tional wave function must be generalized to reflect the 
anisotropy of the effective mass. To take this anisotropy 
into account, we will construct the trial function in the 
form similar to the hydrogenic function of Eq. (f^Oj) . 
but with the appropriately scaled exponential and pre- 
cxponential factors. Due to the axial symmetry of the 
Hamiltonian ([5]) at £ = the £-field vector can be cho- 
sen to lie in the (xz)-plane without any loss of generality. 
Thus our final variational function reads: 

F +z {r)=F z {£) (1 + Q±x + Q\\z) exp(- e /a ± ), (32) 

where Q± = q±± + q 2 ±g, Q\\ = qi\\ + q 2 \\Q , g = 

^Jr 2 + {a\/a 2 — l)z 2 , and F z (£) is the normalization 

constant. Similar variational functions with normaliza- 
tion constants Fj(£) can be defined for any valley sj. 

Friesen investigated the Stark effect for a phospho- 
rous donors in silicon 2 ^ using a single-valley varia- 
tional method and a degenerate perturbation theory (see 
Eq. (|38j) below) to account for the valley-orbit effects. 
Our trial function ([32]) will reduce to that of Friesen's 
if we force parameters q 2 to zero, i.e. set and fix 
q 2 ± = q 2 \\ = 0. As we have seen from Eq. pip it may 
lead to some inaccuracy in the low-filed limit. In what 
follows we will use our improved single- valley variational 
function (|3"2"1) and adopt Friesen's treatment of the valley- 
orbit effects. First, we compute the expectation energy 
of the single valley states, neglecting central cell contri- 
butions. For given 7 and £ we take the expectation value 
of energy E S j to be a function of the parameters a^, and 
the various q\ v and q 2r/ with T) = {||, _L}. The energy 

E sj (a v ,q ln ,q 2 ,; % £) = (33) 

\-fsj I rsj) 

is minimized with respect to these parameters: 

0S(7,£) = dE(n,E) = dEfr,£) = Q 
da v dq lr , dq 2v 

This condition implicitly defines the variational parame- 
ters as functions of the electric field. Using this knowl- 
edge, we can expand the single- valley ground state energy 
in Taylor series around zero field. 



A simple system to test the variational function (|32j) , is 
a hypothetical hydrogenic donor with the isotropic effec- 
tive mass, m± = m\\, i.e. 7 = 1. Since the unperturbed 
Hamiltonian displays spherical symmetry, we can take 
the z-axis to lie along the field direction. If we define our 
atomic units of energy and distance as Eq = h 2 /2to_l<Zq 
and clq = h 2 K/2m±e 2 , respectively, and minimize the en- 
ergy according to (|34p we obtain in the low-field limit: 
Qi± = <72± = 0, g^n = e£/2E , and q 2 \\ = e£/Aa E Q , in 
other words, we recover the exact expansions (|29|) - (p0| . 
where we replace as — > &o and E Ry — > Eq. Similarly for 
the Friesen case with q 2 \\ = we obtain a less accurate 
expansion (|3 1[) . which does not comply with Rayleigh- 
Schrodinger perturbation theory. 

It is expected that our variational function will re- 
turn better corrections than previous methods, especially 
when used for small fields. For silicon, the effective mass 
anisotropy parameter 7 = 0.209. We can no longer freely 
rotate the coordinate system due to a fixed heavy mass 
axis, and we will be required to keep track of the compo- 
nents of electric field parallel and perpendicular to this 
axis. Consider an envelope function of the valley sj, 
F S j(r,£). To study the valley-orbit corrections due to 
the central cell contact potential we will need the values 
of F s j(0,£) = Fj(£) at r = 0. Our computations pro- 
ceed as in the isotropic hydrogenic case, and we obtain 
the second-order expansion of the energy and the central 
cell amplitudes: 

E = Eo-± X \\tf\-\x±£±, (35a) 
FjW/Fo = 1 - \f^£l - \f { l ] £l. (35b) 

Here, £y and £± are the field components parallel and 
perpendicular to the heavy- mass axis of the valley sj, 
Fq = Fj(0) is the normalization constant at zero field, 
and the numerical values of the susceptibilities and coef- 
ficients J*- 2 -* are: 

X|| = 1.71 i ucV(kV/cm)- 2 , (36a) 

X_l = 3.63^cV(kV/cm)- 2 , (36b) 

/,[ 2) = 1.53 x 10~ 4 (kV/cm)" 2 , (36c) 

/} 2) = 2.82 x 10~ 4 (kV/cnVr 2 . (36d) 

We see good agreement with the susceptibilities pre- 
sented in Eqs. (|25a[) and (|25b[) . Additionally, we can 
check the answers obtained using Friesen's previous re- 
sults. As before, by setting q 2 ± and q 2 \\ to zero and hold- 
ing them fixed, we obtain equivalent results for Friesen's 
variational function: 

X\\ = 1.58AieV(kV/cm)- 2 , (37a) 

X± = 3.17 i ucV(kV/cm)- 2 , (37b) 

/,[ 2) = 1.54 x 10^ 4 (kV/cm)- 2 , (37c) 

/{ 2) = 2.69 x 10~ 4 (kV/cmr 2 , (37d) 
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FIG. 3. (Color online.) The plot on the top shows the 
comparison between our results (solid lines) and Friesen's 
(dashed) for the spectrum of a phosphorous donor in silicon. 
The second plot shows the spectrum of a lithium donor, where 
the dashed lines correspond to Friesen's variational function 
(<?2X = ga|| = 0) for Li. 

With the variational functions at hand we can con- 
struct the valley-orbit Hamiltonia n 20 i 37 in the presence 
of the electric field: 

H vo {£) = J2 i E 0i + A 0l ) \si) (si\ + J2 A u \ si ) (- s *\ 

s,4 i,s 

+ A 2y -(l-^)| sl )( s 'j|, (38) 

where E^ may be written as £?oi = Eq — i%j_£ 2 + 
\ (x_l — X||) + ••• ■ Other corrections to the terms in- 
volving Eoi will be left to the next section. The matrix 
elements A q of the central cell contact potential between 
the valley-orbitals © describe the central cell shift (Ao) 
and the valley-orbit splittings (Ai and A 2 ) of the energy 



levels. Following Friesen's prescription^ we parametrize 
the matrix elements A q as follows: 

A 0J = u Ff(£), (39a) 
A y = viFf{£), (39b) 
A 2l] =v 2 E l {£)F J {£). (39c) 

The parametrization in Eqs. (|39a[) - (|39c[) assumes that the 
central cell potential has a contact (i.e. (5-function like) 
form. The fitting parameters depend on the choice of 
the variational functions to ensure that the experimen- 
tal spectrum at zero field is reproduced correctly. For 
our variational functions defined by Eq. (|32[) . it gives 
the following values of the parameters vi in Si:Li: vq=- 
34.54 eV-A 3 , ^=7.09 eV-A 3 , and u 2 =7.09 eVA 3 . In 
Fig. |31 we compare the spectra of P and Li shallow donors 
calculated with ours and Friesen's variational functions 
for an external electric field in 001 direction. For phos- 
phorous impurity our approach does not produce sub- 
stantial differences at low and intermediate fields. For 
lithium, on the other hand, the effect is more significant 
due to the nearly degenerate ground state. 

At this juncture, we need to consider the two com- 
peting effects caused by the electric field. The first ef- 
fect is the direct quadratic Stark shift of the single- valley 
Coulomb binding energies . This effect will be similar 
to that of strain upon the donor spectrum, as we discuss 
below in section lTVl The second effect is based on the fact 
that by pulling the lithium donor electron away from the 
central cell, the electric field reduces the magnitude of 
the matrix elements A q . This will bring the levels closer 
to their "center of gravity" and narrow the overall energy 
spectrum of the Is manifold^ For Li donors subject to 
relatively low electric fields below ionization threshold, 
this spectrum narrowing effect may be treated separately 
and independently from the Stark shift of the valley en- 
ergies Eoi by means of the second-order expansion of the 
variational parameters, as detailed in Section IVl and Ap- 
pendix [B] 

IV. EFFECT OF STRAIN AND ELECTRIC 
FIELD 

We wish to consider the interplay of the Stark effect 
with the effects due to other influences such as strain 
and magnetic field. Both the electric field and strain in 
the silicon lattice will cause a change of the six single- 
valley energies that emerge on the diagonal of the valley- 
orbit Hamiltonian (|3"5]l . By virtue of the ellipsoidal valley 
symmetry these energies will not depend on index s and 
for a given valley si, we can write both the quadratic 
Stark effect and the strain corrections to the energy as 
E sl = E 0i = E - 5i(£) - Si(e jk ) such that 

Si(S) = -\x^ 2 + \ (XI -X||) £?, (40a) 
Si(ejk) = Zd(e X x + e yy + e zz ) + E u eu, (40b) 
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where ejk are the components of the strain tensor at the 
donor location and S<j and S u are the dilation and shear 
deformation-potential constants of Si conduction band 
minimumjSS The nearly degenerate ground state mani- 
fold of the Li donor is very sensitive to strain. However 
a simple dilation merely shifts the levels by the same 
amount, AE = (E d + E u /3)(e xx + e yy + e zz ), and docs 
not alter their separation. To change the relative posi- 
tions of the energy levels and lift the fivefold degeneracy 
of the ground state, at least one of the two linear combi- 
nations of the strain tensor components 

1 \/3 
&9 &zz ~^(&xx ~t~ ^"yy)i ^£ cy (&xx ^yy) (^-0 

must be different from zero^ 

In a similarity, a uniform electric field can also lift the 
degeneracy. The quadratic Stark effect can be described 
by the variables 

w e = 1 -(3nt-l), uj e =^(nt-n 2 y ), (42) 

where n x<v ^ z are components of the unit vector n in the 
direction of the electric field 

£ = £n, n = (n Xl n y ,n z ). (43) 

Therefor to account for the strain effects, we must add 
the strain energies (|40b[) to the diagonal of the valley- 
orbit matrix (see Eq. (|38p) in \si) basis. To the second 
order in £ , we may combine the strain and Stark terms 
and separate them from the zero field valley-orbit Hamil- 
tonian ©. Fixing the average energy of the Is manifold 
at zero, the combined Hamiltonian, describing strain and 
electric field effects, may be written as 

H s = Zu (v 9 Vg + v e fc) , (44) 

where E u = 11.4 eV. Here Vg and V e are operators in the 
space spanned by the six symmetrized orbitals 

Vg =i (\E e )(Eg\ - \E e )(E e \ - \T 2x ){T 2x \ - \T 2y ){T 2y \) 

+ \\T 2z )(T 2z \ + ^ {\A 1 ){E 9 \ + miA^) , (45a) 

V e =-1= (\T 2x )(T 2x \ - \T 2v )(T 2y \) 

+ i (V2\Ai){E e \ - \E e ){E t \+H.c) , (45b) 

where H.C. denotes hermitian conjugation. 

In Eq. (|44[) . the variables vg and v e are linear combina- 
tions of terms describing the effects of strain and electric 
field: 

vg = eg + n£ 2 wg : (46a) 
u £ = e e + k£ 2 u> £ , (46b) 
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FIG. 4. (Color online.) Dependence of the 5 lowest orbital 
energy levels of Si:Li donor electron on the uniaxial [001] 
strain. Dominant characters of the states are given with let- 
ters. Color of the lines also encodes the dominant characters 
of the states for each value of eg . 

where 

and x = X± ~ X\\ = 1-9 (UcV^kV/cm) -2 is the anisotropic 
part of the quadratic Stark effect susceptibility. The 
quantities v e and vg in Eqs. (|46[) can be viewed as "ef- 
fective strain" variables. The quadratic Stark effect due 
to the field lkV/cm is equivalent to a very small strain 
~ 10~ 7 . Larger electric fields of the order of 3kV/cm 
will be equivalent to the strain ~ 10 -6 . 

Elaborating on the analogy between strain and electric 
field, the two effective strain parameters determine the 
orbital response when actual strain or electric fields are 
present. For simplicity, let us initially consider the case 
when v t = but vg is nonzero. As follows from Eqs. (|41 [) . 
(|4*2"|) . (|46a[) and (|46b[) . the degeneracy of the Li ground 
state quintet can be partially lifted either with uniax- 
ial strain along [100], or with biaxial isotropic strain in 
the (xy) -plane, or with some combination thereof. The 
states T 2x , T 2y and E e , which do not contain any contri- 
bution from I iiz), valley-orbitals will remain degenerate. 
The nature of the strain determines the structure of the 
ground state. A tensile uniaxial strain along [001] has 
vg > and favors the (T 2x , T 2y , E e ) triplet ground state. 
Contrariwise, the non-degenerate ground state T 2z and 
the first excited state Eg correspond to a compressive 
uniaxial [001] strain with vg < 0. These findings are 
summarized in Fig. 2J 

One can relate the effective strains produced by an im- 
posed external stress using the three-dimensional Hooke's 
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law for an isotropic material: eg = (1 + v){o zz — \{o xx + 

o yy ))jE and e e = &(1 + v)(o xx - cr yy )/E. Here, £7 
and v are the Young's modulus and Poisson's ratio, re- 
spectively, and o jj is the applied stress component along 
the jth axis. If we wish to create a condition with 
eg 7^ and e e = 0, we require a xx = a yy and find 
eg = (1 + v){p zz — a xx )/E. In short, to prepare a three- 
fold degenerate ground state (T 2x , T^j,, i? e ), a uniaxial 
tension along [001] is sufficient. As we will see below, an 
additional stress along [100] or [010] will result in v e ^ 
and completely lift the orbital degeneracy. 

We can create similar conditions with the electric field. 
To keep v c = at nonzero field, the electric field must 
be confined to the plane formed by the [001] axis and 
either the [110] or the [110] axes. The uniaxial strain can 
be replicated by electric fields lying in this plane. An 
electric fields parallel to [001] replicates tensile uniaxial 
strain, while fields parallel to either [110] or [110] repli- 
cate compressive uniaxial strain. Pointing the electric 
field away from these axes, but still within the plane in 
question, only results in reducing the effective strain v e , 
which vanishes entirely for a field parallel to [111]. 

One can expect a non-trivial interplay of the Zeeman 
and Stark effects due to multiple level crossings in the 
ground state manifold. The possibility of the ground 
state Stark splitting makes the Li impurity unique among 
other shallow donors in Si and opens exciting opportu- 
nities for electrical manipulation of the spin qubits. As 
we will show below, the values of the electric field in 1-3 
kV/cm range are sufficient to produce large changes in 
^-factors near the points of avoided crossing between the 
donor electron energy levels controlled by Zeeman and 
spin-orbit interaction. 

To examine this interplay further, we want to introduce 
a second source of effective strain v e , in this case such 
that \vg\ » \v e \ > 0. A nonzero strain v e separates the 
triplet manifold (T2 X , T2 y , E e ), allowing us to maximize 
the spin-orbit effects and control g-factors. A biaxial 
strain anisotropy or an electric field lying in the (xy)- 
plane will produce nonzero v e . In general, we should 
note that the increase of strain or electric field along a 
single crystallographic axis will change the value of vg by 
v e V3/3. For our purposes, we consider vg w 10~ 4 , which 
will not be seriously influenced by v e ~ 10 -6 . As long 
as \v e \ << \vg\, the correction to vg will not introduce 
any serious errors to the spin-orbit interactions within 
the ground state triplet {T 2x , T 2y , E £ ). 

As before, we can think of the electric field as anal- 
ogous to strain. In general, electric fields lying in the 
(xy)-planc will produce a nonzero effective strain v e . The 
influence of the electric field on v e is maximal along the 
crystallographic axes and vanishes entirely along [110] or 
[110]. The field along [100] produces "tensile" effective 
strain v c > while the field along [010] produces a "com- 
pressive" effective strain v c < 0. 

Despite having similar qualitative features, the electric 
field effects are considerably weaker then those of strain. 
Additionally, strong electric fields will ionize the Li donor 



and may even induce electrical breakdown within the sil- 
icon. For this reason we will consider only moderate elec- 
tric fields in 1-3 kV/cm range and envision their role as a 
means of precise fine tuning and control of the Li donor 
spectrum. We must rely on tensile strain eg to split the 
Li quintet and isolate the (T 2x , T 2y , E e ) ground state 
triplet. The fine tuning by the electric field will then 
induce much smaller effective strain v e . By inducing the 
splitting of the (T 2x , T 2y , E e ) triplet, the electric field can 
be used to explore the effects of the spin-orbit interaction 
through g-factor control of ESR spectra. 



V. EFFECT OF SPECTRUM NARROWING 

Electric fields can influence the donor spectrum in a 
strain-like way, however, higher-order confounding ef- 
fects emerge through the electric field dependence of the 
valley-orbit matrix elements Aoj, Ay, and A 2 ij- Gen- 
erally speaking, these spectrum narrowing effects violate 
relationship between strain and electric field. However, 
the analogy described in Sec. IIV1 remains intact if wc 
are only concerned with the influence of the electric field 
upon the isolated {T 2x , T 2y , E e ) manifold. 

AppcndixlBldetails a procedure for constructing the or- 
bital Hamiltonian with strain, electric field and spectral 
narrowing effects taken. Let us consider a truncated ver- 
sion of the Hamiltonian given in Eq. (|B8[) restricted to the 
{\E C ) , \T 2x ) , \T 2y )} Hilbcrt subspace. We assume that a 
uniaxial tensile stress has separated the other levels and 
isolated the ground state triplet manifold. Then, as de- 
scribed in appendix [B] the spectrum narrowing Hamilto- 
nian will have the form 

Hen = (qe + q v ) (\E e ) (E e \ + \T 2x ) (T 2x \ + \T 2y ) (T 2y \) 



^(\T 2x )(T 2x \-\T 2y ) (T 2y \) 



The form of each q will be 



2F*(v - vi) (/i 2) - /J 



£ 3 W. 



(48) 

(49a) 
(49b) 



% = ^oV (/f + /,i 2) + (/f - /J 2) ) »2) ^ (49c) 

where the parameters i>k, f^ 2 ' and Fq are defined in Sec- 
tion imi 

Because the ground state triplet manifold is well sepa- 
rated from the rest of the states with higher energy their 
contribution to the triplet level shifts is negligible. As 
previously, we eliminate the terms shifting the triplet lev- 
els by an equal amount and set the energy origin at the 
"center of gravity" of the triplet manifold. Then the pri- 
mary effect of spectrum narrowing will consist in renor- 
malization of the Stark susceptibility. At the same time, 
the analogy between electric field and strain remains in- 
tact within the ground state (i.e. triplet) manifold. The 
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new effective strain variable controlling the splitting of 
the triplet levels can be written as 



1^2 
V e = e e + -— O W e 



(50) 



where \' is the effective susceptibility given by 

X ' = x + 2^(,o-- 1 )(/f-/ | i 2) ) 

« 1.5/ieV(kV/cm) _2 . (51) 

While this susceptibility is slightly reduced, we demon- 
strate that it is still large enough to induce dramatic 
shifts of electron g-factors as well as overall reshaping of 
ESR lines. 



VI. EFFECTS OF ZEEMAN AND SPIN-ORBIT 
INTERACTIONS 



The full donor electron Hamiltonian, 

H = H V0 (£,e ij ) + H z + H s , 



(52) 



includes the valley-orbit, strain and Stark effects as well 
as the spin-Zeeman and spin-orbit interactions^ char- 
acterized by (/-factor anisotropy and by two spin-orbit 
constants^ 

g± = 1.9984, = 1.9994, (53) 
Ai = 2.6 ficV, A 2 = 6.9 ^cV. (54) 

We can recast the spin-orbit Hamiltonian using vector 
operators similar to those introduced by Watkins and 
Ham£: 



H so — -(Ai-fci + A2-L2) • 



where 



Li = - -= \ si ) s l n i x n j] s '( s 'j\ 



(55) 



(56) 



ij ss' 



~* ij ss' 

(57) 

where r = (1,1,1), <x = (a x ,& y ,& z ), and &i are the con- 
ventional Pauli matrices. Similarly, the Zeeman Hamil- 
tonian can be expressed as: 



er B + eV] \is)&iBi(is\ 



(58) 



where e = (g\\ - g±)/g± = 5 x 10 4 . 

In ESR studies, pulses of ac magnetic field excite donor 
electron spins coupled via the magnetic dipole interaction 
to a cavity mode with fixed frequency ujq. By varying 



the strength of the external static magnetic field B, the 
cavity mode is excited every time when the resonance 
condition, E n — E m = Huiq, is fulfilled. Here E n — E m 
is a Zeeman splitting for the transition between states 
with predominantly opposite spin orientations. The spin- 
orbit interaction^ gives rise to shifts hA nm of the Zeeman 
splitting for spin-flip transitions 



E„ - E m = g±/J, B B + hA, 



(59) 



Each transition corresponds to a specific value of static 
magnetic field B nm at which the resonance with cavity 
mode is achieved. It is customary to formally introduce 
^-factors instead of B nm for each transition 



9n 



HbB, u 



E n — E n 



(60) 



where /is is Bohr magneton. In general, the Zeeman 
shift A nm is a function of magnetic field magnitude (B) 
and orientation (r) as well as of the value of the effective 
strain (v e ) at the donor location. Introducing a cyclotron 
frequency, u> = j-g±HBB, the values of g-factors g nm at 
each resonance equal 



9n 



g±- 



UJQ 



(61) 



where lu = uj nm is a root of algebraic equation (cf. (|59[0 



u = lu - A nm (w,T, v e ). 



(62) 



It is of interest to study the g-factor for each transition as 
function of magnetic field orientation and effective strain: 



1 



flu 



91- 



i(w»m,T,!) e ) 



(63) 



Assuming that the bare Zeeman energy is much greater 
than the spin-orbit interaction, hu^> Ai,A2, the corre- 
sponding shifts in Zeeman splitting A nm (uj) for each 
transition can be found from the successive orders of per- 
turbation theory corrections to the energy levels E n . As 
a result, each shift A nm can be obtained in a form of an 
expansion in inverse powers of u: 

A nm = d™ m + dr^- 1 + d% m uj- 2 + d^ m oj- 3 + ... . 

(64) 

Because the energy corrections contain the powers of bare 
Zeeman splittings in the denominator, the expansion (|64[) 
does not have positive powers of w. The term with zero 
power of lu in (|64| occurs when states with the same 
spin orientation have very close energies and coupled by 
spin-orbit interaction. The latter then splits the Zeeman 
transition frequencies already in zcroth order in Ai ; 2 (this 
happens in our case as will see below). Plugging (|B"4")) 
into (|62[) and inverting the series expansion, one can find 
a root of (|62p for a given transition, tjj=U) nm , in terms 
of inverse powers of tu . Then, expressing magnetic field 
B nm in (|60[) via io nm gives the ^-factor in the form 



9 = 91 



dn 

1 + — 

LOQ 



dl +di dl + Modi + d 2 



(65) 
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where we omitted the state indexes for simplicity. In 
general, for a given Li donor the g-factors for individ- 
ual transitions depend on both the electric field and the 
strain at the donor location through the effective strain 
variables vg, v e , as well as on the magnetic field orienta- 
tion. 

Typically, <?-factor shifts are very small in the ls(Ai) 
ground state of substitutional shallow donors in silicon 
because the spin-orbit interaction only involves excited 
states ls(E + T 2 ) and the corresponding corrections con- 
tain very large energy denominators for the valley-orbit 
splitting between the singlet and the rest of Is levels. In 
the Li donor, the situation is different because the bundle 
of levels, ls(E + T 2 ), with little or no valley-orbit split- 
ting is now a ground state. It is exactly the zevo-th order 
term oc do in Eq. (|64p. resulting from this degeneracy, 
that gives rise to very large g-factor shifts. Below we will 
analyze various types of the spin-flip transitions in the 
vicinity of the avoided crossing. At the avoided crossing 
the situation is quite involved because the spin-orbit cou- 
pling mixes the states with different orbital characters. 
The role of (xy)-plane strain or electric field is to lift the 
level degeneracy and to promote direct spin-flip transi- 
tions between the states with the same orbital characters. 



VII. G-FACTOR CONTROL WITH ELECTRIC 
FIELDS 

Here we consider a uniaxial tensile stress along [001] 
amended by an electric field along [100]. This combina- 
tion of stress and electric field leads to effective strains 
vg, v e such that \vg\ >> \v e \. The small, yet nonzero, 
effective strain v e splits the triplet levels according to 
Eq. (|45b[) . If we include the electron spin in the triplet, 
we have a sextet of spin-orbital electron states. For uni- 
axial tensile stress 0001 ~ 30 MPa and magnetic field 
0.343 T (corresponding to an ESR cavity mode of fre- 
quency ujq = 9600 MHz), the energy gap separating 
a ground state spin-orbit sextet from the higher-lying 
states is of the order of 1 meV (cf. Fig. [4]), and it is 
much greater than the Zeeman splittings within the sex- 
tet (~ 40 /i.eV). Therefore, for sub-Kelvin temperatures, 
the system can be described by a reduced Hamiltonian 
operating on the sextet of electron spin states. In what 
follows, we shall study the effects of spin-orbit interac- 
tion, small i> £ , and magnetic fields upon this sextet. 

The dependence of the lower triplet of energy levels 
on the effective strain variable v e is shown in Fig. [5] for 
the case of a magnetic field in [001] direction. The level 
structure is symmetric with respect to ±i> £ . The behavior 
of these levels can be understood from the fact that when 
stress and magnetic field are aligned with the same crys- 
tal axis (z) the eigenstates of the donor electron Hamil- 
tonian split into two subspaces. The first subspace is 
formed by the states \E € , i), |T 2a ,,— \T 2y ,—\) while 
the other one corresponds to the opposite spin projec- 
tions on magnetic field direction. The states belong- 
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FIG. 5. (Color online.) Plot shows the level diagram in the 
vicinity of one of the avoided crossing at v e —0. Blue, green 
and red line colors correspond to the dominant orbital char- 
acters T2 y , Ec, and T2 X , respectively. Symbols indicate 
the eigenstates (|70[) for the corresponding energy levels. 
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FIG. 6. (Color online.) Plot shows dependence of the six 
lowest energy levels on [100] electric field with uniaxial tension 
0001 = 30 MPa and magnetic field B001 = 0.343 T. 

ing to different subspaces are eigenstates of the operator 
Z = R z (ir)a z , which commutes with the Hamiltonian. 
Here R z (tt) is the operator of rotation through angle ir 
about the z-axis and Z\<&±) = ±|$±). Due to this sym- 
metry, the eigenvalues of Z are good quantum numbers 
and the spin-orbit interaction couples the states within 
each subspace, but not between the subspaces. 

For the eigenstates of the donor Hamiltonian, neither 
the orbital characters (T 2x , T 2y , E e ) nor the spin projec- 
tions (f, 4-) are good quantum numbers. Nonetheless, we 
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will label these eigenstates as T 2x ^, E e ^ etc., keeping in 
mind that this is predominant, albeit approximate, char- 
acter of a given eigenstate. 

As shown in Fig. [31 when a [100] electric field increases 
from zero, the levels T 2x and T 2y arc shifted quadratically 
with the field (linearly with v e ) in opposite directions 
while the level E e does not change in the first order in 
v t (its variation occurs in the second order from coupling 
to higher lying states Eg, A\ as AE = —n t v1 with r\ t = 
45. keV). For a given magnetic field and effective strain 
v e = Ve such that 



where the parameter 



(66) 



(ve « 6.3 x 10" 6 for B 00 i = 0.343 T), there exists an 
avoided crossing between the levels corresponding to the 
states with dominant characters l^i/j^) and \E e , — 
Higher in energy by approximately .g/is-B there exists 
another avoided crossing between \E e , ^) and \T2 X , — \) 
(£ ~ 11 kV/cm in Fig. [6]). For v t ~ — w e (w), there exists a 
similar pair of avoided crossings where the orbital charac- 
ter T 2x is replaced by T 2y and vice versa. We note that for 
the twice smaller value ofv e «3xl0~ 6 (£ ~ 8 kV/cm), 
the energy levels of the states \T 2x , — |), |T 2y , ^) undergo 
a real crossing because these states are not coupled to 
each other by the spin-orbit interaction. 

A different type of avoided-crossing exists for v e ~ 0. 
It occurs between the levels corresponding to the states 
with the same spin orientation and predominant orbital 
characters \T 2x ) and \T 2y ). We note the energies of the 
states T 2x , T 2y , and E e with the same spin projection 
are very close to each other. Despite this, the state E e 
is not coupled to the other two by the spin-orbit interac- 
tion. Therefore energy levels of the states with dominant 
characters | , ± -i ) are only weakly perturbed by spin- 
orbit interaction and effective strain v e (see above). 

The level splitting for the the avoided crossing between 
states with opposite spin orientations, A2, is about twice 
larger than that between the states with the same spin 
orientation (corresponding to Ai). However, in either 
case the, level repulsion at the avoided crossings is much 
smaller than the Zeeman energy, and therefore, the anal- 
ysis of the avoided crossings can be done to a leading 
order within the two-state approximation. In what fol- 
lows, we will use the two-state approximation and con- 
sider the effects of the spin-orbit interactions near the 
avoided crossings at v e = and v e ~ 6.3 x 10~ 6 . 

We consider truncated Hamiltonians to describe the 
avoided-crossings at v e —Q- As E c is not coupled by 
the spin-orbit interaction to the other states, a two-level 
Hamiltonian in the basis of the corresponding pairs of 
states \T 2xi \), \T 2y ,\) and \T 2x ,-\), \T 2y ,-\) is ap- 
propriate. The form of these Hamiltonians can be written 
as 



Hid 

' 2 = ± T 



1 






W ±4 \ 











(67) 



2E U 

W = —V t 



Aiv/3 



(68) 



controls the detuning from the avoided-crossing reso- 
nance due to strain and/or Stark effects. The energies of 
the pairs of states near the two avoided crossing points 
equal 



huj ctAi 



(69) 



where cr=±l is a new orbital quantum number (in ad- 
dition to the spin z-projection ±5). The corresponding 
cigenfunctions are 

|* ff ,±i) = (±i C(T \T 2x ) + Vl - 4 |T 2j/ )) 1±|> , (70) 

where 

r i V 2 

a w 

V2 I \/l + w 



±1. 



(71) 



It is of interest to consider the magnctic-dipolc transi- 
tions between states with opposite spin orientations near 
the avoided crossing. The matrix elements of the Pauli- 
matrix a x for spin-flip transitions between states with 
the same orbital number a equal 



a w 



Vl + w 2 ' 



±1. 



(72) 



It is seen from that the frequencies of these transi- 
tions (as well as that of the transition between the states 
\E e , ±1/2)) are very close to w forming a triplet of center- 
lines of magnetic dipole transitions. The matrix elements 
for spin-flip transitions between states with opposite or- 
bital numbers a = ±1 equal 



1 



(73) 



According to (|69j) , the frequencies of these transitions are 
offset from w by ±7i _1 Ai\/l + w 2 , forming a doublet of 
satellite lines. 

Away from the avoided crossing, \w\ 3> 1, the ma- 
trix elements for transitions between states of the same 
a dominate. This behavior can be understood from the 
fact that, in those regions, each of the hybridized or- 
bitals ([70)l is dominated by a single orbital character (cf. 
Fig.E^b)). Specifically, when — w S> 1, 



\*-,±±)^\T 2y ,±±), |*+,±§) 
and when 



|T 2x ,±i>, (74a) 



|*-,±i)^|T 2 .,±i), |*+,±I)^|T 2 ,,±I). (74b) 

Therefore, away from the avoided crossing (i.e. for 
\w\ 3> 1) the spin-flip transitions (|72[) conserving a are 
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FIG. 7. Matrix elements for the spin-flip transitions vs. 
electric field near the avoided crossing at u 6 =0. Dashed lines 
(almost undistinguishable in the plots) show the analytical 
results (|72[) . (I73[) . The lines marked by: (1) correspond to 
matrix elements for the two direct transitions; (2) correspond 
to the matrix element for the transition within the Zeeman 
doublet \E t , ±1/2} which is very close to 1; and (3) correspond 



to the matrix elements for the two satellite lines ^ 



t 



between states with the same dominant orbital character, 
and the corresponding matrix elements approach unity. 
At the same time, the transitions in (|73[) corresponding 
to satellite lines, arc between states with increasingly or- 
thogonal orbital characters as w increases and therefore 
the corresponding matrix elements approach zero in this 
limit. 

The behavior of the matrix elements for spin-flip tran- 
sitions at the avoided crossing (w = 0) is very different 
than that away from it. It follows from (fT0| that at w=0 
the pairs of states |*+, +|) and l* - ,-^), 
|^E' + ,+i) are, respectively, symmetric and antisymmet- 
ric superpositions of the orbital states i\T2 X ) and \T2 V ): 



■2y> 



V2 

-i\T 2x ) + \T 2 y) 



V2 



l±*>. 
*l±§> 



(75a) 
(75b) 



The spin-flip transitions (|72[) that are dominant away 
form avoided crossing \w\ 3> 1, connecting the same 
orbital states, become suppressed at w = because 
they connect the symmetric and antisymmetric super- 
positions of i\T 2x ) and \T 2y ). At the same time, the spin- 
flip transitions (|73[) connecting the states with different 
orbital characters, which are suppressed away from the 
avoided crossing {\w\ 3> 1), become dominant at w = 
where they connect the superpositions of i\T 2x ) and \T 2y ) 
with the same symmetry. This behavior is evident from 
Fig. (J7J giving the dependence of matrix elements on elec- 



tric field (or effective strain v e ) obtained using the exact 
numerical solution for the eigenstates of the donor elec- 
tron Hamiltonian. 

There exist four distinct spin-flip transitions between 
the states of a spin-down doublet ^ ■ and those of a spin- 
up doublet . We shall denote the corresponding g- 
factors as g a ^ a i (a, cr'=±l). Based on Eqs. (see also 

Fig. O there are two satellites corresponding to #-,+, 
g~,+ and two closely spaced center lines corresponding 
to .9-.-, g+.+- There is a third center line with g-factor 
gE e corresponding to the transition between the states 
|-E e ,±i). All g- factors can be obtained by numerical 
solution of Eq. (|62|) for the corresponding transitions. 
The results are shown in Figs. [5J 

The numerical results can be very closely approxi- 
mated analytically when the (xy)-plane strain splitting of 
T 2x , T 2y is much smaller than the Zeeman splitting, that 
is, for |« e |<g;« e (jSrJ)) . We find the coefficients in by 
a perturbation theory expansion in the spin-orbit inter- 
action constants Ai,2 using the basis of "correct" states 
([70]) in zeroth order. All five g-factors will have the form 



g = g Q (l + S) + A*, 

1 f^2_ 

8 \hu!o 



3 V 5o 



(76a) 
(76b) 



where go is a singlet g-factor and Al corresponds to dif- 
ferent shifts for the various ^-factors. 

For the satellite lines, the shift A^ at equals up to the 
3 rd order in X 2 /ujq: 



g XiD(v e ) 5 \i\ 2 2 D(v e ) 

' ~\~ (7 ■ 

hujQ — a\iD(v e ) 8 (huo) 3 



(77) 



Here, we used a secular approximation in (|63[) keeping 
the terms linear in Xi/lu in the denominator because 
they arc proportional to the dimcnsionlcss parameter D, 
which increases away from the avoided crossings 



D(V C ) = VT+ 



w- 




(78) 



where w is given in (|68|) . At the avoided crossing, ^-factor 
shifts from the center line g${l + 5) for satellites are near 
2X 1 /fuv ps 0.13 (for uj = 9600 MHz). As we discussed 
previously, while the p-factor shifts for satellite lines in- 
crease away from the avoided crossing, the strength of the 
lines decrease. Note, however, that with detuning from 
the avoided crossing corresponding to £ — 3.5kV/cm 
(v e ps 10 -6 ) the matrix element is only suppressed by a 
factor of 4 while the g-factor shift is already of the order 
of unity. This giant change in ^-factor can be seen by 
comparing the curves <?-,+ and <?+.- in Fig. |8£b), and 
the curves labeled 3 in Fig. 

Fig-UJa) shows the p-factors <?+,+, <7-,- and gE e corre- 
sponding to three very closely spaced center lines. Their 
splitting is not resolved within the 2-level picture near 
the avoided crossing given by Eqs. (|69p . Wc must take 
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FIG. 8. (Color online.) g-factors vs electric field near the 
avoided crossing at v e —0. Figure (a) shows three weakly split 

center lines, g++,g and gE t - Centerline go(l + 5) (I76al ) 

is shown in bold. Figure (b) shows the satellite lines g±,^. 
Electric field 3.5 kV/cm corresponds approximately to the 
strain v e « 10~ 6 . Dashed lines (unresolved in some plots) 
correspond to analytical approximations given in Eqs. (|76a|l . 
(77), (79) and (80). 



into account the transitions between the states 1^^) and 
the corresponding states \E e ) with opposite spin ordina- 
tion, to determine the splitting. The values of A c give 
the shifts of three centerlines from the common center. 
For g ajy (cr=±l) they are obtained from the expansions 
of A£. in powers of u>q ■ 



\l{D-a){D*\l + {D-cT)\l) 



Ah 2 u$D 



(79) 



Here a = ± and D (78) is a dimcnsionlcss detuning from 
resonance. Note a slight asymmetry between the line 
shifts that occurs in the 4th order in X2/HUJ0. Analytical 
expressions match very closely exact numerical results as 
demonstrated in Fig. (5Ja) where analytical results are 
shown with dashed lines. The g-factor shift for the tran- 



sitions between the \E e ,±l/2) states equals 



AK 2 lo 



2X1) 



8h 4 u 



(80) 



Note that the dependence of this ^-factor on the effec- 
tive strain variable v e occurs only in the 4 th order in 
\2/hiuo leading to its slight decrease with v e also shown 
in Fig.lSJa). In the leading order, g-factor shifts are sym- 
metric: ~ — Ai. At the avoided crossing (D = 1), 
^-factor shifts Al and A^, nearly coincide as can also be 
seen in Fig.[5]Ja). Away from the avoided crossing D^> 1, 
the lines A± asymptotically approach each other as they 
are shifted down by {\ 2 / tvjOo) 2 / A m 0.0076 from the line 
Ae c - Wc finally note that in the expressions above (77), ( 
[75), (SO), we have replaced the singlet g- factor g Q with 
the integer value 2. 

The above predictions can be better understood and 
visualized if we plot the ESR lineshapes for different elec- 
tric fields. Microwave fields stimulate the spin-flip transi- 
tions, with resonance emerging as a decrease of the inten- 
sity of the received microwave field. The ESR lineshape 
results are predicted through a simple Lorcntzian model. 
To account for thermal excitation, wc assume that the to- 
tal contribution to the ESR lineshape from a single donor 
exposed to microwave light of frequency luq will take the 
form of 



x{B 1 £ 1 e)=Z- 1 Y, 



e -E q /k B T 



qr 



(B,£,e), (81) 



q,r>q 



where Z = ex P(~ E n /kBT). Here the line- 

shape of a single spin-flip transition between the states 
denoted as q and r. The single-transition contribution 
reads 



x qr (B, £, e) = 



1 



2n r 2 



A r ,(B,£,e)] 



(82) 



where A rq = (E r — E q )/H and Xqr = (q\ o'x \ r ) is the 
spin-flip transition amplitude between these two states. 
The parameter T expresses the degree of line-broadening 
due to dispersive mechanisms such as thermal effects and 
spontaneous phonon emission. 

Away from the circled avoided crossings in Fig. [B) the 
transitions arc strongest between the states of the same 
orbital character. At an electric field of 3kV/cm, we 
would expect to see three strong absorption lines cor- 
responding to the three resonant transitions between 
|T 2x ,t) o |T 2x ,i), |T 2y ,t) o \T 2y ,i), and \E 6 ,t) 
\E e ,l). At the avoided crossing with v e = 0, a very 
different situation will be in effect. The transition cor- 
responding to Tix and T 2y is suppressed, and new tran- 
sitions with greatly shifted g-factors will emerge. In the 
above nomenclature, the transitions corresponding to the 
(/-factors g+,+ and <?-,- are suppressed in favor of satel- 
lite transitions corresponding to (/-factors <?+,- and <?-,+■ 
The satellite transitions are very sensitive to the random 
strains, which implies that the likelihood of their detec- 
tion is extremely low. For this reason we will concentrate 
on the centerlines and their electric field dependencies. 
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FIG. 9. (Color online.) Modeled ESR lineshape (9600 MHz 
cavity mode, 1.3 gauss FWHM, 4 Kelvin) for different £-fields 
along [100]: (a) £=0, (b) £= 1.5 kV/cm, and £=3 kV/cm. 
The emergence and shifts of and T-2 y ESR lines are clearly 
seen. 



The g-factor control emerges as the ability to selec- 
tively turn these lines "on" and "off", and to shift the 
lines themselves. Fig. [TJ] shows the predictions of our 
model. The suppression of the two centerlines T 2y and 
T 2x at zero field, i.e. near the avoided crossing v e = 0, 
and their emergence at higher fields can be clearly seen. 
The zero-field suppression of the lines is due to the spin- 
orbit interaction. The two lines emerge because the elec- 
tric field shifts the energy levels away from the avoided 



crossing where the influence of the spin-orbit coupling is 
much weaker. The overall dependence of the intensity of 
T 2x and T 2y lines on the electric field closely follows that 
of the spin-flip matrix element shown in Fig. [7] As we 
see in Fig. [9j in the absence of random strain the Stark 
effect induces a dramatic shift of the ESR lines on the 
order of 10 gauss. 

In Figs [10] and [HI we show the effects of random strains 
upon the ESR lineshapes and predicted Stark shifts of Li- 
donor ^-factors. We assume that the experimental ESR 
signal >c can be modeled as an ensemble average over the 
Gaussian distribution of the random strains: 



ji(B,£) = J n(e)>c(B,£,e) de. 



(83) 



Here n(e) is the Gaussian distribution function and the 
integration is taken over the strain variables eg and e e . 
We further assume that a strong uniaxial tensile stress is 
applied to the sample, and the random internal strains 
shifting eg will have a negligible contribution to the over- 
all ESR signal, i.e. may be safely ignored. Thus we 
consider only the effect of random variations of e e , as- 
suming they are described by the Gaussian distributions 
with different standard deviations (uncertainties) Ae e . 

To better understand the random strain effects pre- 
sented in Figs. [TOl and [TT1 let us recall our previous find- 
ing that the electric field can strongly affect the Zee- 
man splittings of the centerlines only in the vicinity of 
the avoided crossings (Fig. [5]). In the domain of random 
strains with Ae e not exceeding 10 -7 , most of the donors 
reside near the avoided crossing v e = 0. Thus, for the 
majority of these donors, the centerline transitions T 2y 
and T 2x will be suppressed at £=0; however these transi- 
tions will emerge for £-fields exceeding 1 kV/cm. At £=3 
kV/cm, the donor spectra will be shifted from v e = to 
v e , exceeding the standard deviation Ae e . As a result, we 
see pronounced Stark shifts of the ESR lines on the order 
of 10 gauss (see Fig. [9] and the top panel of Fig. [TO"]) . 

The strain disorder randomly shifts the spectra of indi- 
vidual donors away from the the avoided crossing v e — 0. 
For broad distributions most of the donors are subject 
to the random strains e e that exceed v e corresponding to 
£ =3 kV/cm. The ^-factors of these donors are saturated 
and are insensitive to the applied f- field (see Fig. [5] (a)). 
As a result, the majority of the Li spins do not display 
any significant line shifts induced by the electric field, 
and the Stark features are washed out of the average sig- 
nal of the ensemble. As we see from the bottom panel 
of Fig. [H for Ae e = 5 ■ 10" 7 , the two lines T 2x and T 2y 
already exist at £ =0, which means they are induced by 
the random strain. The latter also broadens the lines 
and pins them down at the zero-field positions. In the 
domain of larger random strains Ae e > 5 • 10~ 7 the two 
lines will merge, broaden and eventually collapse, result- 
ing in a spectrum insensitive to the electric field, with one 
narrow line E e . The electric field dependencies of the av- 
erage ^-factors corresponding to different Ae e (Fig. [TTj) 
clearly display flattening with increase of Ae e . 



17 




3390 3400 3410 3420 3430 3440 
Magnetic Field (gauss) 



280 



c 200 
o 

'I 

o 

< 120 



40 









T 2 y 












I 


1 


Tlx 













3390 3400 3410 3420 3430 3440 
Magnetic Field (gauss) 

FIG. 10. (Color online.) ESR absorption lineshapes near 
v e = for different levels of the random strains and different 
electric fields £ . The solid and dashed lines represent the 
lineshapes at £ = and £ = 3 kV/cm respectively. The 
Stark shifts are prominent for Ae s = 10~ 7 (top panel) but 
are almost entirely washed out for Ae £ = 5 • 10 -7 (bottom 
panel) . 



The typical values of the random strain uncertainties 
Ae e , which allow for observation of appreciable Stark 
shifts in Li-doped Si, should be in the range of Ae e ~ 
10~ 7 (Fig. fTUj) . While this requirement is very stringent, 
the possibility of the growth of low-random-strain Si ma- 
terials has been demonstrated in photoluminescence ex- 
periments^ with 28 Si epilayers grown on natural silicon 
substrates. Yang et al£2- were able to resolve splitting of 
bound exciton lines due to the lattice constant mismatch 
Aa/a ~ 10 -6 between the epilayer and the substrate. 
The observed exciton linewidths are at least an order of 
magnitude smaller than the splitting, which is indicative 
of the required level of the random strains ~ 10 -7 . The 
latter are most likely caused by isoelectronic impurities or 
complexes (e.g. carbon^), which implies that the chemi- 
cal purity of the material is a key factor in reducing these 
strains. 
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FIG. 11. (Color online.) Stark shifts of electron g-factors 
for different levels of the random strain. The solid lines are 
the exact p-factors at zero random strain, shown before in 
Fig. [8] The dot-dashed lines correspond to Ae E = 10~ 7 and 
the dashed lines represent the case of Ae e = 5 • 10 -7 . 



VIII. SUMMARY AND DISCUSSION 

We have developed a theory for the Stark effect for 
lithium donor spins in silicon. The anisotropy of the 
effective mass leads to the anisotropy of the quadratic 
Stark susceptibility, which we determined using the 
Dalgarno-Lewis exact summation method^. The the- 
ory is asymptotically exact in the field domain below Li- 
donor ionization threshold, relevant to the Stark-tuning 
ESR experiments^. Using this theory as a calibration 
tool we devised a new variational wave function for a 
shallow donor in the electric field. The variational func- 
tion replicates the exact small-field asymptotic results 
and is robust for large fields. 

With the calculated Stark susceptibilities and the new 
variational function at hand, we predicted and analyzed 
several important physical effects. First, we observed 
that the energy level shifts due to the quadratic Stark ef- 
fect are equivalent to and can be mapped onto those pro- 
duced by an external stress^. Second, we demonstrated 
that the Stark effect anisotropy, combined with unique 
valley-orbit splitting of a Li donor in Si, spin-orbit inter- 
action and specially tuned external stress, may lead to 
a very strong modulation of the donor spin ^-factor by 
the electric field. Third, we investigated the influence of 
random strains on the g-factor shifts and quantified the 
random strain limits which are necessary to observe the 
ESR-Stark shifts experimentally. 

The ability to control ^-factors with electric field 
and/or stress is a crucial component of many solid-state 
based quantum computing schemesi^— Of particular in- 
terest for QIP proposals is the situation that emerges 
when states with opposite spin-orientations are allowed 
to mix through the spin-orbit interaction leading to the 
avoided crossings. This has implications not only for g- 
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factor control but also for long-range inter-donor clastic- 
dipole coupling, suggesting the possibility of control- 
lable interactions between isolated donor qubits.— Both 
of these capabilities can be utilized to implement a uni- 
versal set of gates based on the Ising Hamiltonian4i This 
would require ability to grow Si multilayer structures 
with a given value of the in-plane strain. The promis- 
ing technique, which can be used for practical purposes 
of achieving precise control of the uniaxial and/or biaxial 
strain, is the growth of Si films on compliant surfaces.— 
Another interesting possibility is the utilization of piczo- 
actuators to implement local stress control of individual 

I 

We seek to solve the differential problem 



impurities 
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Appendix A: Derivation of Eqs. (|19[1 



dl 



/-A 



dz 



f 



where £ = (x,y, za»/a±). First, we use the product rules 



V 2 (cos(<t>)f(r,6)) =cos<f> ( V 



r 2 sin 2 i 



2m±a 3 L 



f(r,0), 



dz ' 



(Al) 

(A2a) 
(A2b) 



and similar rules for the derivatives of sin0 • f(r,9). If we substitute / in the form of Eqs. (|T5)) into Eq. (|A1|) and 
use the rules (|A2|) we are able to separate the (f> dependence and obtain partial differential equations for /y (r, 0) and 
f±(r,8) as follows: 



1 



2 • 2 

r z sin 



V 2 - 2 
2 



o_ 

Or 
d_ 

Or 



f\\- X 
/±-A 



a 2 
dz' 
if 



z d 

2 ~7T 

r dz 

z d 



/ll 
f± 



-rsinO. 



(A3a) 
(A3b) 



^dz 2 r dz / 

These equations can be further simplified by employing product rules involving the associated Legendre polynomials 

1(1 + 1) 



r 2 sin 9 



/(r)P z m (cos 0) = P[ n (cos0) 



and substituting Eqs. ([17]) into Eqs. (|A3|) . This yields: 



Y,Pi(cos0) (r 2 D r - 1(1 + 1)) fy = 

i 

2 P, 1 (cob^) (r 2 D r - 1(1 + 1)) f ±il = 



Xr 2 53 D z (/|| ti Pi (cos 0))-r 3 cos 
i 

Xr 2 53 D z (f ±tl Pf(cos0)) + r 3 sm0 



(A4) 

(A5a) 
(A5b) 



where D r and D z arc defined in Eqs. ([T5 



If A = 0, Eqs. (|A5j) have particular solutions /|| = 
f(r) cos0 and f± = f(r) s'm(0). Indeed, we obtain an 
ordinary differential equation for f(r) 



I 

which has the particular solution 



(A7) 



^+»(l-r)|-V 



(A6) 



The operator D z mixes the Legendre polynomials with 
different I and does not allow to separate r and 0. How- 
ever, the problem can be reduced to a chain of coupled 
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ordinary differential equations. To proceed, we employ 
another set of rules: 



f) flf f f)P m 

-f(r)Pr(cos9) = cos^P/"^ - L sinW^p 



cos(9)P l 



I + m 



PP. 



i 2[ + 1 i-i 



(A8a) 



9P ; ™_ (Z + l)(Z + m) i(i-m + l) 



m 



21 + 1 



2Z + 1 

(A8c) 



We use these rules to rewrite r 2 D z /;(r)P; m (cos 9) as 
r 2 D z (/,(r)Jr(cosfl)) =fe a d£ a + 



+ ^ 2 7J? 2 )/Kr), (A9) 



where a™, /3; m , 7 ; m are second order radial differential op- 
erators having a common structure of Eq. (|20j) with the 



coefficients af^ , and 7^ given explicitly as 



(0) 



(Z + m+ 1)(Z + m + 2) 
(2Z + 3)(2Z + 5) 



aS = -2(Z + 3) 
= + + 3) 



(0) 



/3 (1 



2Z 2 + 2Z - 2m 2 - 
(2/-l)(2/ + 3) 



l.m 



4m 2 



1 



2/ 2 + 2Z - 2m 2 



(0) _ (j - - m - 1) 
7 '< m ~ (2Z — 1)(2Z — 3) 

7S=2(/-2) 
7£ = «*-2) 



We then multiply both sides of Eqs. (|A5|) by 
P z m (cos 6<) sin 6> and integrate d6»sin(6')P/ n (cos6 ) )(...), 
to eliminate ^-dependence at the expense of mixing // 
with different values of I. This procedure leads to 
Eqs ([n 



(AlOa) 

(AlOb) 
(AlOc) 
(AlOd) 

(Alia) 
(Allb) 
(Allc) 
(Alld) 



(A12a) 

(A12b) 
(Al2c) 
(A12d) 



Appendix B: Spectrum Narrowing Effect 

The matrix of the central cell potential in the basis of 
valley-orbitals reads: 



H' 



1 A 0x 


Ala 


A-2xy 


A-2xy 


&2xz 


A 2xz 


Aix 


Aoa; 


A-2xy 


^2xy 


&2xz 


A 2a;z 


A2xy 


A 2x y 


Aoy 


Ai„ 


A 2 yz 


A 2yz 


&2xy 


A^xy 


Ai, 


Ao, 


&2yz 


A 2yz 


&2xz 


&2xz 


A 2 yz 


A 2j(2 


A 0z 


A u 


\&2xz 


&2xz 


A 2yz 


A 2 y z 




A 0z 



(Bl) 



Transforming this the symmetrized-orbital basis, we find 



(D A 


.4 


D 








o\ 


A 


D E e 


c 











B 


C 


D E e 




















D, 




















Dy 





V 

















(B2) 



Here the energy zero is shifted to the "center of gravity" 
of the manifold. Explicit expressions for the diagonal 
matrix elements are 



D A =- (A lx + Ai y + A u + 4A 2xy 



4A 



2xz 



D E 6 = 



~ 6 (A - 



Dec = 



Ao y — 2Aq z — Ai x — Ai y 
4A 2x y + 8A 2 xz + 8A2y Z ), 

A 0y - 2A 0z + 3A lx + 3A ly - 



6 (A °* 



D x — — (2Aoa; — Aqj, — Aoz — 3Ai x ) , 
Dy = q (2Ao a — Ao x — Aoz — 3Ai y ) , 



v 
D- 



(2A 0z - A 



O.r 



A 



0y 



3A lz ). 



H 4A 2yz ) , 
(B3a) 

4A lz 
(B3b) 

12A 2a j / ) , 
(B3c) 

(B3d) 
(B3e) 
(B3f) 



Similarly, the non-diagonal matrix elements can be ex- 
pressed as 



A = - ^ (A 0x 



A 0y - 2A 0z + A lx + A ly - 2A lz 
AA 2xy - 2A 2xz - 2A yz ) , (B3g) 

Ala; _ Aly 



B =-^j- (A 0x — Ao y 



C = -^{A Qx 
6 



A lx - Aly 



2A 2x2 — 2A 2j/z ) , 

(B3h) 

-4A 2x z +4A 2yz ). 

(B3i) 



At zero electric field, F x = F y = F x = Fq. Using these 
values in Eqs. (|39ap ~ p9cp . and substituting them into 
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the matrix (|B2|) , we find the zero- field valley-orbit Hamil- 
tonian to be 

















\ 





D (o) 

E 




















D (o) 

E 




















4°) 




















4 0) 





^ o 














4°V 



(B4) 



where each term on the diagonal is such that 

Df =F*(v x + 4v 2 ), (B5a) 
4 0) =F*(v l -2v 2 ), (B5b) 
4 0) = -Fgvx. (B5c) 

We identify each term on the diagonal with the singlet, 
doublet, and triplet binding energies of the Is donor man- 
ifold. Since these values are well known from experiment, 
we can use the constants vo, V\, v 2 to reproduce the spec- 
trum of any shallow donor in Si. 

The spectrum narrowing Hamiltonian is defined by 
H sn = H vo — . To evaluate the extent of spectrum 
narrowing due to displacement of the electron away from 
the donor site, we write 

F j =F (l + 6 j ), (B6) 

where Sj is a function of the field. The values of Sj can 
be obtained variationally. Our analysis yields 

^ = ^(/f -/ifVl—/?^ (B7) 

where £j is the component of the field lying along the 
ith valley. For fields up to lOkV/cm, the values of 5 



remain below 0.02. In the valley-orbit matrix the cou- 
plings are A oj = v F§(l + 2Sj), A tj = ^F 2 (l + 25j), 
A-2ij — v 2 F§ (I + Si + Sj) where the terms to the first 
order in 5 (second order in £ ) are retained. 

To simplify the matrix expressions, we express for our 
valley-orbit matrix in terms of Si. Then the spectrum 
narrowing Hamiltonian matrix can be written as the sum 
of three simpler matrices 

H sn = F 2 (s e Ug + 5 c U t + <y/„) . (B8) 

The parameters 5^ , determining the strength of the nar- 
rowing, are given by 

S g =1 (-S x -S y + 2S Z ) , (B9a) 

Se ~ (S x - S y ) , (B9b) 

S n =S x + S y . (B9c) 

The component matrices may be expressed in the in the 
symmctrized-orbital basis 

U e =(>i + 4^ 2 ) + (vq + 2i/! - 4z/ 2 ) \E ) (E e \ 

+ (2^ - 3n) \T 2z ) (T 2z \ + v a \A X ) (E e \ + H.C. 

- v (\E C ) (E e \ + \T 2x ) (T 2x \ + \T 2y ) (T 2y \) (BlOa) 

U e =V3 {ya - vi) (\T 2x ) (T 2x \ - \T 2y ) (T 2y \) 

- v e \E ) (E e \ + v a \Ai) (E e \ + H.C. (BlOb) 

U v = -vi (\T 2x ) (T 2x \ + \T 2y ) (T 2y \ + \T 2z ) (T 2z \) 
+ (v 1 -2v 2 )(\E e ){E e \ + \E e )(EA) 
+ (v 1 +Av 2 )\A 1 )(A 1 \ (BlOc) 

with v a = y2 {yo + v\ + v 2 ) 1 and vg = Vq + V\ — 2v 2 . 
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